Iterative real-time path integral approach to nonequilibrium quantum transport 
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We have developed a numerical approach to compute real-time path integral expressions for 
quantum transport problems out of equilibrium. The scheme is based on a deterministic iterative 
summation of the path integral (ISPI) for the generating function of the nonequilibrium current. 
Self-energies due to the leads, being non-local in time, are fully taken into account within a finite 
memory time, thereby including non-Markovian effects, and numerical results are extrapolated both 
to vanishing (Trotter) time discretization and to infinite memory time. This extrapolation scheme 
converges except at very low temperatures, and the results are then numerically exact. The method 
is applied to nonequilibrium transport through an Anderson dot. 

PACS numbers: 73.23.-b, 73.23.Hk, 72.10.-d, 02.70.-c 
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I. INTRODUCTION 

Quantum transport has attracted theoretical and ex- 
perimental research since it offers the possibility to inves- 
tigate quantum many-body properties at as well as out 
of thermodynamic equilibrium [l| . The ongoing improve- 
ment in miniaturization down to the nanometer scale al- 
lows to study electron transport in ultra-small devices, 
e.g., in single molecules or artificially designed quantum 
dots @, H, 0, H, @, 0| ■ There is a broad variety of interest- 
ing physical effects, due to interactions or the nonequi- 
librium conditions arising when a bias voltage is applied 
to the source and drain electrodes [1, 0, [lj| • These range 
from Coulomb blockade via coherent (e.g., resonant tun- 
neling) transport to the Kondo effect, to name but a few. 
While on the experimental side, progress stems from an 
increased control of fabrication processes, many theoret- 
ical works deal with refined approximation schemes ap- 
plicable in different parameter regimes. However, exact 
theoretical results — either analytical or numerical - 
for nonequilibrium quantum transport systems are rare, 
mainly because of the lack of adequate methods allowing 
to tackle such questions. There clearly is a considerable 
need for numerically exact methods to describe nonequi- 
librium quantum transport, both to check analytical (and 
usually approximate) approaches and to connect theory 
to experiment. Here, we mean by nonequilibrium trans- 
port specifically those phenomena which go beyond the 
standard approach of linear response to the applied bias 
voltage. 

In this work, we propose a novel numerical scheme de- 
noted as iterative summation of real-time path integrals 
(ISPI), in order to address quantum transport problems 
out of equilibrium. Many-body systems driven out of 
equilibrium are known [ll|, GJ, [Of to acquire a steady 
state that may be quite different in character from their 
ground-state properties. Details of the steady state may 
depend on the nature of the correlations, as well as on 
the way in which the system is driven out of equilibrium. 
While there are a variety of nonperturbative techniques 
in place to study equilibrium systems, many of these 
methods cannot be extended to nonequilibrium systems 



in a straightforward way. 

Different approximations have been pursued previously 
in order to tackle nonequilibrium situations. For in- 
stance, transport through an Anderson dot in the Kondo 
regime has been described theoretically along several dif- 
ferent lines, e.g., for the asymptotic low-energy regim e by 
Fermi liquid theory [l4[, via interpolative schemes [l5| . 
using intcgrability concepts |16l|. or by the perturbative 
renormalization group [12L Il7| . Transport features of the 
Anderson model have also been discussed by perturba- 
tion theory in the interaction strength [lH, [3] . Sophis- 
ticated nonperturbative methods have been developed 
in order to extract exact results out of equilibrium for 
special models, where intcgrability is available. For in- 
stance, the interacting resonant level model presently en- 
joys much interest [20, SH, [H, H, SI SI Q . Univer- 
sal aspects of nonequilibrium currents in a quantum dot 
have been discussed in Ref. [2l|. In various perturbative 
regimes, this model has been addressed by field theory 
methods [H, SI SE SI] . However, the underlying the- 
oretical concepts are still under much debate, and par- 
tially conflic ting results were reported in the literature, 
cf. Refs. HIHIH!. 

The above discussion shows that precise numerical sim- 
ulation tools arc in great demand in this field. Let us 
briefly review the existing numerical methods available 
for nonequilibrium transport. Much effort has been di- 
rected to the application of renormalization group (RG) 
techniques to this problem. For instance, a perturba- 
tive real-time RG analysis of nonequilibrium transport 
has been performed [27|, and steps towards the appro- 
priate generalization of the density matrix renormaliza- 
tion group technique have appeared pel l29l | . In addition, 
the functional RG approach has been generalized to non- 
equilibrium recently [30| , It forms a systematic and fast 
perturbative expansion scheme, yielding reliable results 
for small interaction strengths, where the infinite hierar- 
chy of equations for the self-energy can be cut after the 
first few steps. Moreover, a possible extension of Wil- 
son's numerical RG approach to nonequilibrium has been 
discussed [3l|, SH] . Another line of attack considered nu- 
merically exact real-time quantum Monte Carlo (QMC) 
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simulations (for the corresponding equilibrium case, see 
Refs. [lii HH). Due to the sign problem, however, these 
calculations become increasingly difficult at low temper- 
atures. Refined multi-level blocking techniques [HI, HH 
allow to achieve further progress, but numerical simu- 
lations remain hard within this approach. (For recent 
progress, however, see Ref (3tJ •) Let us also mention the 
flow equation method, which has been applied to study 
the nonequilibrium Kondo effect [H| . Finally, a very re- 
cent development considers non-standard ensembles to 
describe steady state transport, but involves a numeri- 
cally troublesome analytic continuation [3!| . 

Our ISPI approach, described in detail below, provides 
a novel alternative, and in principle numerically exact, 
method to tackle out-of-cquilibrium transport through 
correlated quantum dots. The method is deterministic 
and, hence, there is no sign problem. It is based on 
the evaluation of the full nonequilibrium Keldysh gen- 
erating function, including suitable source terms to gen- 
erate the observables of interest. It builds on the fact 
that the fermionic leads induce self-energies that are non- 
local in time, but which decay exponentially in the long- 
time limit at any finite temperature. Thus, a memory 
time exists such that within this time span, the corre- 
lations are exactly taken into account, while for larger 
times, the correlations can be dropped. This allows to 
construct an iterative scheme to evaluate the generating 
function. An appropriate extrapolation procedure allows 
then to eliminate both the Trotter time discretization er- 
ror (the Hubbard-Stratonovich (HS) transformation be- 
low requires to discretizc time), and the finite memory- 
time error, yielding finally the desired numerically exact 
value for the observables of interest. The extrapolation 
procedure is convergent for not too low temperatures, 
since then memory effects are exponentially small for 
long times. If the extrapolation converges, we thereby 
obtain numerically exact results for nonequilibrium quan- 
tum transport properties of interacting systems. 

The ISPI scheme is implemented here for the example 
of the well-known single-level Anderson impurity model 
[40L Ull . [H, [H, HH, but appropriately modified, it can 
be applied to other quantum dot models as well. The 
nonequilibrium current is calculated from a generating 
function, represented as a real-time path integral in the 
Keldysh formalism. After a Hubbard-Stratonovich trans- 
formation employing an auxiliary Ising spin field, all 
fermionic degrees of freedom (of the dot and the leads) 
can be integrated out exactly, however, at the price of 
introducing time-nonlocal self-energies for the leads and 
a path summation over all Ising spin configurations. For 
this, an iterative summation scheme is constructed, ex- 
ploiting that the time-nonlocal correlations in the lead 
self-energy effectively decay exponentially at finite tem- 
perature T or bias voltage V, thereby setting the char- 
acteristic memory time. For larger times, the correla- 
tions arc exponentially small and will be neglected. The 
full time-discretized Keldysh Green's function (GF) of 
the dot then assumes a band matrix structure, where 



the determinant can be calculated iteratively using its 
Schur form. The scheme is constructed such that within 
the range set by the memory time, the path integral is 
evaluated exactly. The remaining systematic errors are 
the Trotter time discretization and the finite-memory er- 
ror. Both, however, can be eliminated in a systematic 
way based on a Hirsch-Fye type extrapolation proce- 
dure. Based on the above construction principle, our ap- 
proach cannot be applied at very low energies (T, V — > 0), 
where, fortunately, other methods are available. For fi- 
nite temperatures, the requirement of not too long mem- 
ory times can be met, and the spin path summation re- 
mains tractable. 

The present paper is organized as follows. In Sec. [U 
we introduce the model for the quantum dot coupled to 
normal leads. We present the computation of the gener- 
ating function for the nonequilibrium Anderson model in 
Sec. [Hi] in the presence of an external source term, which 
allows to calculate the current as functional derivative. In 
Scc. lIVI wc introduce the numerical iterative path-integral 
summation method, from which we obtain observables of 
interest. We give a detailed discussion of the convergence 
properties of our method and describe the extrapolation 
scheme in Sec. [V] For several sets of parameters, we will 
present results in Sec. EH followed by a discussion in Sec. 
IVIII Throughout the paper, we set h = fee = 1. 

II. MODEL 

We consider the Anderson model [I(| given by the 
Hamiltonian 

T~t = Hdot + H leads + H-T 

= ^2 Eoafio- + Un^hi + ^(efep - Vp)c\ p<J c kpa 

- E [*p4p«A + H • (!) 

kper 

Here, E 0a = E + oB with a = j, j= ± is the energy of 
a single electron with spin a on the isolated dot, which 
can be varied by tuning a back gate voltage or a Zee- 
man magnetic field term oc B. The latter is assumed 
not to affect the electron dispersion in the leads. The 
corresponding dot electron annihilation/creation opera- 
tor is da-/d a , with h„ = d\d a with eigenvalues n a = 0, 1, 
and U denotes the on-dot interaction. For later pur- 
pose, it is convenient to use the operator identity h-\hi = 
\{n^ +ni) — \ {h] — hi) 2 , thereby introducing the shifted 
single-particle energies e^ a = Eo a + U/2, which yields the 
equivalent dot Hamiltonian 

Hdot = H dot fi +Hjj = e 0(T n<T - — (n-|- - hi) 2 . (2) 

In Eq. fT]), Ck p denotes the energies of the noninteract- 
ing electrons (operators Ckpa) in lead p = L/R = ±, 
with chemical potential /i p = peV/2. Dot and leads 
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are connected by the tunnel couplings t p . The observ- 
able of interest is the (symmetrized) tunneling current 
I=(Il- Ir)/2, 



I(t) 



kpa 



(3) 



where I p (t) = -eN p (t) with N p {t) = (£ ka c\ pcT c kpa ) t . 
The stationary steady-state dc current follows as the 
asymptotic long-time limit, / = lim t ^ 00 /(t). We have 
explicitly confirmed that current conservation, Il + Ir = 
0, is numerically fulfilled for the ISPI scheme. In the 
presence of a finite bias voltage, F / 0, the Keldysh 
technique [45|, [46|, [47| provides a way to study nonequi- 
librium transport. In this formalism, the time axis is 
extended to a contour with a = ± branches, see Fig. [IJ 
along with an effective doubling of fields. The Keldysh 
GF is Gf{t ai t' p ) = -i{Tc[^i{t a )^){t' p )]), where T c de- 
notes the contour ordering of times along the Keldysh 
contour, and i,j = L, R, correspond to fields represent- 
ing lead or dot fermions, respectively. We omit the spin 
indices here, remembering that each entry still is a di- 
agonal 2x2 matrix in spin space. In general, the four 
Keldysh components are linearly dependent, such that 
G ++ + G~~ = G H + G h . However, the commonly 
used Keldysh rotation [i^, [47j seems to offer no advan- 
tages here, and is not employed in what follows. 



III. GENERATING FUNCTION 

Let us then discuss the generating function, which con- 
tains all relevant information about the physics of the 
system. First, we want to integrate out the lead fermion 
fields, and, in addition, perform a discrete Hubbard- 
Stratonovich transformation. This allows to integrate 
out the dot fields as well, and we are then left with a 
discrete path summation. We start with the fermionic 
path-integral representation of the generating function, 



m-J 



V 



~\_d a ,da-, 



Ckpa i &kpa 



(4) 

with Grassmann fields (d, d, c, c) (Note that we use the 
same symbols for fermion operators and Grassmann fields 
throughout the paper for better readability). We intro- 
duce an external source term S^ (defined below in Eq. 
©), which allows to compute the current at measure- 
ment time t m . 



I(t m ) 



-«J^hiZ[r/] 



(5) 



r/=0 



We note in passing that it is also possible to evaluate 
other observablcs, e.g., the zero- frequency noise, by intro- 
ducing appropriate source terms and performing the cor- 
responding derivatives. The action is S — Sdot + Si ea d s + 



St + S v , where 
Sdot = Sdot.o + Sjj 
dt 



U 



Sieads 

Sj • 



^ d a (idt - too) d a + ^-(rJ| - nj.) 2 



/ dty^c kpa (id t - e kp + [i p )c kpa , 
Jc 



kpa 



<c 
ierj 



p^kpa^a H - h.C.j 



^ p{ip^kpad(j ^pdaC-kpa ) (^m ) • 
kpa 



(6) 



The interaction term Su in Eq. ([6|) does not allow to di- 
rectly perform the functional integration. Apart from 
this, all other terms are quadratic in the Grassmann 
fields and thus define the 'noninteracting' (quadratic) 
part. Note that the level shift +U/2 is here included 
in the noninteracting sector. 



A. Noninteracting part 

Before turning to the interacting problem, we briefly 
discuss the Keldysh GF of the noninteracting problem. 
Let us first integrate out the lead degrees of freedom. We 
remain with an effective action for the dot, which is non- 
local in time due to the presence of the leads. In partic- 
ular, after Gaussian integration, the generating function 
reads 



Zni[rf\ 



V 



~[da,d„ 



ASdot.O + Senv) 



(7) 



with 



Senv = f dt f dt'J2da(t)l^ L (t,t')+^R(t,t') 
J C J (J „ \ 



+ ! fW,*')-7H(i,*')] 

X [S(t-t m )-S(t' -t m )]\da(t') 



(8) 



For the source term, the physical measurement time t m 
can be taken at one branch of the contour, see Fig. [l] As 

we fix t m on the upper (+) branch, the ( ) Keldysh 

element of the source term self-energy vanishes, see the 
time-discretized version below. The "f p in Eq. ([8]) repre- 
sent the leads, and their Fourier transforms are explicitly 
given as 2 x 2 Keldysh matrices, 



1p 



M = i?v ( o /( " 9 7/ p) 7 \ 9 7? /( " -f>\ ) ■ 0) 

p \ 2 - 2f(w - fi p ) 2f[u - n P ) - 1 J 



Here, we have used the fact that the leads are in thermal 
equilibrium, f(uj) = l/(e"/ T + 1). Moreover, we take 
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the standard wide-band limit, with a constant density of 
states per spin channel around the Fermi energy, p(ep), 
yielding the hybridization T p = irp(eF)\t p \ 2 of the dot 
level with lead p. For the sake of clarity, we assume from 
now on symmetric contacts, Yl = Yr = Y/2. The gen- 
eralization to asymmetric contacts is straightforward. In 
the next step, still for U = 0, we may also integrate over 
the dot degrees of freedom. We obtain the nonintcracting 
generating function 

ZniM = II dct [-iG^(t,t') + r)i; J (t,t')] , (10) 



The solution is A± = A' ± iX" with 



5 t X' = sinh -1 y/sm{6 t U/2), 5 t X" = sin" 1 y/sin (6 t U/2). 

(13) 

Note that the overall sign of X± is chosen arbitrarily, 
but does not influence the physical result. Unique- 
ness of this Hubbard-Stratonovich transformation re- 
quires USt < 7T. To ensure sufficiently small time dis- 
cretizations, in all calculations one should then obey the 
condition max(U,e\V\, \e \,T) < l/5 t . 



where G ^(t,t') follows from 

G Qa (uj) = [(oj - e 0a )T z - 7l(w) - 7ij(w) 
1 



r2 + ( W - e 0CT ) 2 

w - e 0a + iY(F - 1) iYF 

iT(F-2) -u> + e 0a + iY(F - 1) 

where t z is the standard Pauli matrix in Keldysh space, 
and F = f(w + eV/2) + f(u> - eV/2). Moreover, the 
self-energy for the source term is obtained as 

£ J (M') = e -[idt,t')- lR (t,t')] 

x [6(t - t m ) - S(t' - t m )] . (11) 

Up to this point, we have discussed the noninteracting 
case. In order to treat interactions, we now perform a 
Hubbard-Stratonovich transformation. 



B. Hubbard-Stratonovich transformation 

Let us therefore turn back to the real-time action Sdot 
of the dot in the presence of interactions, U ^ 0, see Eq. 
©. For later numerical purpose, it is beneficial to pro- 
ceed with the discussion in the time-discretized represen- 
tation of the path integral [48[ . In order to decouple the 
quartic term, we discretize the full time interval, t — NS t , 
with the time increment St- On each time slice, we per- 
form a Trotter breakup of the dot propagator according 

t0 e iS t {H +H T ) = e iS t H T /2 el S t H 0e iS t H T /2 + Q^2^ wheye 

H = Hdot + Hi ea ds- By this, we introduce a Trotter er- 
ror 0, [5(], [HI which, however will be eliminated from 
the results in a systematic way p2\ |53|. see below. Next, 
we use a discrete Hubbard-Stratonovich transformation 
[5tl . l54l , [55j for the interacting part, which introduces 
Ising-like discrete spin fields = (s+ , s~ ) on the a = ± 
branches of the Keldysh contour (with s" = ±1) on the 
71-th Trotter slice. For a given Trotter slice, we now use 



s ± t S t U(n r -n l ) 2 /2 _ 



IE 



D -S t \±s ± (n r -n l ) 



(12) 



s±=± 



For U > 0, noting that n-f — n± = 0, ±1, the sum can be 
carried out and gives the condition 

cosh(<5 t A±) = cos(S t U/2) ±ism(5 t U/2). 



C. Total GF and generating function 

After the Hubbard-Stratonovich transformation, the 
remaining fcrmionic Grassmann variables (da-, d a ) can be 
integrated out at the cost of the path summation over the 
HS Ising spins {s}, 

{«} CT 

with the total Keldysh GF written in timc-discrctizcd 
(1 < k,l < N) form as 



(15) 

where a, j3 = ± labels the Keldysh branches, and the 
noninteracting GF is 



.a/3 



G 



2T 



l5t< - k - l)uJ G 0<7 (Lu). 



(16) 



Note that Go a (t,t') depends only on time differences 
due to time-translational invariance of the noninteract- 
ing part. Moreover, the self-energy kernel stemming from 
the external source term follows as 



J kl 



a/3 _ a/3 
7i M ~ "fR,kl 



[SmkS a ,+ - S m lS/3, + ] , (17) 



where 7 Pi fcz = 7 P (ife — U) and the measurement time is 
t m = mS t . 



IV. THE ITERATIVE PATH-INTEGRAL 
SCHEME 

Let us now exploit the property (see, e.g., Ref. [5^|) 
that each Keldysh component of Goa,M decays exponen- 
tially at long time differences (|A; — l\ — > 00) for finite 
T, see Eq. (fl~6| . We denote the corresponding time 
scale (correlation or memory time) by t c . In Fig. [2] 
we show typical examples of Re G^ kl for different bias 
voltages V. The exponential decrease with time is il- 
lustrated in the inset of Fig. [UJa) where the absolute 
value I Re G^~ kl | is again plotted for the same parame- 
ters, but in a log-linear representation. For large bias 
voltages and low enough temperatures, e.g., at V > Y 
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and T < 0.2T, the decay is superposed by an oscillatory 
behavior, see Fig. [2] Since the lead-induced correlation 
function decays as ~ cos[eV(t — t')/2]/ smh[irT(t — t')], 
the respective correlations decay on a time scale given 
by t~ ~ max^^T 1 , eV). (The correlation function also 
has an additional eo-dependence of the decay character- 
istics.) Thus, the exponential decay suggests to neglect 
lead-induced correlations beyond the correlation time r c . 
This motivates an iterative scheme which exactly takes 
into account the correlations within t c , but neglects them 
outside. The exponential decay has to be contrasted to 
the case T = V = 0, where correlations die out only 
algebraically, and our approach is not applicable. 

Let us then face the remaining path sum in Eq. (fT4|) . 
In the discrete time representation, we denote with to = 
< tjv = NSt the initial and final time, and tk = k5 t , 
see Fig. [1] The discrctized GF and self-energy kernels for 
given spin a are then represented as matrices of dimen- 
sion 27V x 27V ". For explicit calculations, we arrange the 
matrix elements related to Keldysh space (characterized 
by the Pauli matrices r) and to physical times (tj., t{) as 
r ® (k,l). In particular, the ordering of the matrix ele- 
ments from left to right (and from top to bottom) repre- 
sent increasing times. The lead-induced correlations thus 
decrease exponentially with growing distance from the 
diagonal of the matrix. For numerical convenience, we 
evaluate the generating function in the equivalent form 

Z[r,]=^Y,l[detD <7 [{s},r 1 }, (18) 

{s} a 

with D a = Goo-G" 1 . Explicitly, this reads 

(19) 

The normalization prefactor M does not affect physical 
observables, and is put to unity. The time local nature 
of the on-dot interaction is now present only in disguise, 



since the matrix product lets Ising spins appear line- wise. 
By construction, we have to sum over 2N auxiliary Ising 
spins, and the total number of possible spin configura- 
tions is 2 2N . 

Next, we exploit the above-mentioned truncation of 
the GF by setting Dm = for \k — l\St > t c , where 

t c = KS t (20) 

is the memory time, with K the respective number of 
Trotter time slices. The GF matrices then have a K- 
band structure. Note that in the continuum limit, where 
K = N,St — > and N — > oo, the approach is exact. 

To prepare the basis for the iterative scheme, we now 
transform the GF matrix to Schur's form, which then 
allows to calculate the determinant in a straightforward 
way. In general, a quadratic block matrix D given as 

D = ^ ^ d ) ' w ^ a k° m § invcrtible, can be repre- 
sented in its Schur form, i.e., after one step of Gaussian 
elimination. Hence, by multiplying from the left with a 
lower triangular matrix, 

b = LD = ( -m- 1 I m ) ( c d ) = ( d - ca^b ) 

where I n ( m ) represent unit matrices of dimension n = 
dim a (m = dimd), and 6, c do not have to be quadratic 
themselves. Clearly, the determinant is invariant under 
this transformation. The (2, 2) element in this block no- 
tation is often referred to as the Schur complement of 
the matrix D. This representation thus allows to write 
the determinant as det(-D) = det(a) det(d — ca" 1 ^). We 
can now establish the iterative scheme as follows. We 
start from the full GF in the matrix representation as 
defined in Eq. (fT!?)) . After the memory truncation, the 
TV x N— matrix (in time space) assumes a band structure 
represented as 



D 1 



D 1 







D 21 D 22 D 23 



D = D (h 



D 32 D 33 D 34 



D 43 D 44 



£)N K -IN K 



V o 



jjN k N k -1 £)N K N K 
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where the single blocks are K x K— block matrices defined as (I, I' = 1, 



,N 



K 



D 



W 



( D d-i 



(l-l)K+l,(l'-l)K+l 



\ Aif,(i'-i)Js:+i 



D 



(l-l)K+l,l'K 



D 



IK.l'K 



The number N of Trotter slices is always chosen such that Nk = N/K is integer. The elements Dki are given in Eq. 
(|T5|) . with their dependence on the Ising spins s k kept implicit. Each Dm still has a 2 x 2— Keldysh structure, and a 
2x2 spin structure. Then, we rewrite the generating function (|18|) as 

Z[ V ] = J2 deb{D n [*t,...,4]} 

s ± s ± 

x det {d (2iNk) [s± +v . . . , s%] - D 21 [s K+v . . . , sf K ] [D n [sf, . . . , D 12 [s K+1 , . . . , , (21) 

where the Nk — 1 x Nk — 1— matrix D(2,jv r ) is obtained from Dn^^ by removing the first line and the first column. 

In order to set up an iterative scheme, we use the following observation: to be consistent with the truncation of the 
correlations after a memory time KSt, we have to neglect terms that directly couple Ising spins at time differences 
larger than the memory time. This is achieved by setting 



jjl+2,1 + 1 jjl + l.l f£)Ml _1 £)U+1 £)l+l,l+2 _^ q 



(22) 



within the Schur complement in each further iteration step. Note that we do not neglect the full Schur complement 
but only those parts which arc generated in the second-next iteration step. With this, we rewrite the generating 
function as 



N K -1 



£ det 4» I] dct{z?^'+V* + i,- 



1=1 



_ D l+l,U ± ± i 

u l b lK + l> ■ ■ ■ i b (l + l)Kl 



(l-l)K+l> ■ ■ ' ' °IK\ 



D 



i>4 



(23) 



Then, one can exchange the sum and the product, and by reordering the sum over all Ising spins, one obtains 



m = E 



%N K \ S N-K+X ' ' ' * ' S n\ ' 



(24) 



where Zjy K is the last element obtained from the iterative procedure defined by (Z = 1, . 



,N 



K 



Zl+l[ S lK+V ■ ' ■ ' S (i+1)A'] — E ^[ S (i-l)JC+l' ■ ■ ■ ' S IK ' S iif+1> ' ■ ' ' s (Z+l)if]^[ s (;-l)_ff+l' 



*(I-1)K + 1 *IK 



(25) 



The propagating tensor A; can be read off from Eq. (|23[) as 



A, = det{D ,+1 -' +1 [^ +1 , 



_r,Z+Ur ± 



S (Z+l)if] 

(j+i)ifJ p L s (;-i)fsM-i! ■ • ■ i s ik\ 
I 



(26) 



Note that we use here the notion of tensor in the sense 
of a multi-dimensional array and not in the strict mathe- 
matical sense defined by the transformation properties of 
this object, see also Ref. [57j. The iteration starts with 
Z 1 [st,...,s%] = det{D n [4,..., S %]}. 

The current is numerically obtained by evaluating Eq. 
(|5|) for a small but fixed value of rj; we have taken rj = 



0.001 for all results shown below. By this, we obtain the 
full time-dependent current I(t m ) as a function of the 
measurement time t m (0 < t m < NS t ). At short times, 
this shows a transient oscillatory or relaxation behavior, 
which then reaches a plateau value from which we read 
off the steady-state current I. 
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V. CONVERGENCE AND EXTRAPOLATION 
PROCEDURE 

In order to render the scheme exact, we have to elimi- 
nate the two systematic errors which are still present up 
to this point, namely, (i) the Trotter error due to finite 
time discretization St = t/N, and (ii) the memory error 
due to a finite memory time t c = KSf. The scheme 
becomes exact by construction in the limit K —> oo 
and St — » 0. To perform this limit in a straightforward 
way is not possible due to the exponential dependence 
of the array sizes on K. However, unless temperature 
is very low, we can eliminate both errors from the nu- 
merical data in the following systematic way: (i) We 
choose a fixed time discretization St and a memory time 
t c . A reasonable estimate for t c is the minimum of l/|eV| 
and 1/T (see above). With that, we calculate the cur- 
rent I{8t,T c ), and, if desired, the differential conductance 
dl(6t, T c )/dV (the derivative is performed numerically for 
a small AV = 0.01T). The calculation is then repeated 
for different choices of St and r c . (ii) Next, the Trot- 
ter error can be eliminated by exploitin g th e fact that 
it vanishes quadratically for St — ► [49|, l50l . l5l| . For a 
fixed memory time r c , we can thus extrapolate and ob- 
tain dI(r c )/dV = dl(5 t — ► 0, T c )/dV, which still depends 
on the finite memory time t c . The quadratic dependence 
on St is illustrated in Fig.[3]for different values of U. Note 
that each line corresponds to the same fixed memory time 
r c = 0.5/r. (iii) In a last step, we eliminate the mem- 
ory error by extrapolating for l/r c — * 0, and obtain the 
final numerically exact value dl/dV = dl(r c — > oo)/dV . 
For the dependence on l/r c , we empirically find a regular 
and systematic behavior as shown in Fig. [H The r c — > oo 
value is approached with corrections of the order of l/r c , 
see Fig. |4j However, it should be stressed that tempera- 
ture and voltage affect the convergence properties in dif- 
ferent ways, as is already clear from Fig. [2j Importantly 
enough, when T and V are such that the extrapolation 
scheme described above converges, numerical exactness 
is warranted. 

We have implemented the iterative scheme together 
with the convergence procedure on standard Xeon 2GHz 
machines. Computations are then only possible for K < 
7 due to the limited memory (RAM) resources available. 
Typical running times for the shown simulation data are 
approximately 15 hours for K — 5. With a second, paral- 
lelized version of the code, we have been able to take into 
account up to 2 14 summands in Eq. (|25p . corresponding 
to K = 7, within passable running times of a few days. 



VI. RESULTS 

Next, we discuss the results obtained by the applica- 
tion of the iterative procedure to the Anderson model. 
As pointed out before, see Sec. IIII Al we consider the 
symmetric case, Tr = Tl = T/2, with ^l/r = ±eV/2. 
In what follows, we measure energies in units of T. Un- 



less noted otherwise, all error bars for the shown data 
points, which are due to the Trotter and memory extrap- 
olation scheme, are of the order of the symbol sizes in 
the figures. Our scheme yields the full time-dependent 
current I(t), including transients as well as the asymp- 
totic steady-state value. Fig. \5\ shows typical results 
for the current I(t), for two parameter sets, namely (1) 
U = 4,T,eV = 2T,T = OAT (black circles), and (2) 
U = 0.5T,eV = 0.6r,T = T (red squares). The first 
set is of interest in the context of the nonequilibrium 
Kondo effect, where the Kondo temperature (for eo = 0) 
is Tk = y/TU/2exp(— irU/8T); for parameter set (1), 
this yields T K = 0.29r. For eV > T K , analytical re- 
sults for the steady-state current arc available [HI, [58| . 
shown as solid line in Fig. This indicates that the 
ISPI method is capable of approaching the nonequilib- 
rium Kondo problem. For both parameter sets, and for 
many others not shown here, we observe that after a tran- 
sient relaxation behavior, I(t) settles at a plateau value, 
which then defines the stationary current discussed in the 
following. 



A. Validation of the algorithm: comparison with 
exact and perturbative results 

As a simple warm-up check, let us briefly compare our 
numerical results to the exact result for U = 0. Fig- 
ure Ela) shows the stationary current (obtained already 
at t m = 10/r) as a function of eo for T = 0.1T and 
eV = 0.2T. The I-V characteristics is shown in Fig. EJb) 
for T = O.ir and eo = 0. Both cases illustrate that the 
numerical result coincides with the exact one. Other pa- 
rameter sets have been checked as well, and agree with 
the well-known U = analytical solution. 

Second, in order to validate the reliability of our code 
for finite U, we compare the numerical results to a per- 
turbative calculation, where the interaction self-energy 
is computed up to second order in U [H^]. In order to 
respect current conservation, this calculation is possible 
only at the electron-hole symmetric point, eo = B = 
fl8| . First-order terms give then no contribution, and the 
self-energy corresponds to just one diagram [59| ■ For a 
detailed comparison, we plot mostly the interaction cor- 
rections, 5 A = A(U) - A{U = 0), with A being the cur- 
rent J, the linear conductance G, or the nonlinear con- 
ductance dl/dV, respectively. Figure [7| shows the results 
for 51 as a function of the bias voltage for U = 0.1T and 
U = 0.3r. For U = O.ir, we perfectly recover the pertur- 
bative results, which confirms the reliability of our code 
even in the regime of nonlinear transport. Clearly, the 
corrections are small and negative, which can be rational- 
ized in terms of Coulomb blockade physics, as transport 
is suppressed by a finite interaction on the dot. 

For U = 0.3r, the current decreases even more, and 
the deviations between the ISPI and perturbative results 
are also larger. The relative deviation for U = 0.3T is al- 
ready ~ 30 — 35%, illustrating that perturbation theory 
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is already of limited accuracy in this regime. Although it 
well reproduces the overall tendency, there is significant 
quantitative disagreement. The differences arc even more 
pronounced for U = Y, as shown in Fig. [5] for 51 (main) 
and for 8{dl /dV) (inset). Here, second-order perturba- 
tion theory does not even reproduce qualitative features. 



B. Comparison with master equation approach 

Next, we compare our numerical approach with the 
outcome of a standard master equation calculation [4lj . 
The master equation is expected to yield reliable results 
in the incoherent (sequential) tunneling regime, T Y, 
where a description in terms of occupation probabilities 
for the isolated many-body dot states is appropriate. The 
transport dynamics is then described by a rate equation 
for the populations, where the time-dependent rates are 
obtained from lowest-order perturbation theory in Y [4l| . 
The results are shown for U = Y in Fig. [HI both for the 
nonlinear and the linear differential conductance interac- 
tion corrections. As temperature is lowered, interaction 
effects become more important, as seen from the exact 
ISPI results. This corresponds to the emergence of co- 
herence effects for T < Y, which are clearly not captured 
by the master equation in the sequential tunneling ap- 
proximation. However, for T > 4r, interaction correc- 
tions are washed out, and the master equation becomes 
essentially exact, cf. Fig. [51 Similarly, from our ISPI 
results, we observe that interaction corrections are sup- 
pressed by an increasing bias voltage as well. To give 
numbers, for T — 1.25r, we find in the linear regime 
SG = -0.073e 2 //i, whereas S(dI/dV) = -0.062e 2 //i for 
eV = 3Y. 



C. Small bias: eV < V 

For sufficiently small bias voltage, the current is linear 
in V, and we can focus on the linear conductance and its 
interaction correction SG. Figure [TO1 shows SG as a func- 
tion of eq for different magnetic fields B, taking U = Y 
and T = 0.1 (to be specific, we have chosen eV = 0.05r). 
For B = 0, two spin-degenerate transport channels con- 
tribute, and a single resonant-tunneling peak at e = 
results. The interaction corrections are most pronounced 
on resonance, eo = 0. For B ^ 0, the spin-dependent 
channels arc split by Ae = 2B, resulting in a double-peak 
structure. The spin-resolved levels are now positioned at 
eo±B due to the Zeeman splitting. Again the interaction 
corrections are largest on resonance, and the double-peak 
structure of G(eo) is transferred to <5G(eo), cf. Fig. [101 We 
find no evidence for an interaction-induced broadening 
of the resonant-tunneling peak compared to the nonin- 
teracting case. The width of the Lorentzian peak profile 
for B = is determined by Y at sufficiently low T, and 
broadens as T increases. Here, the double-peak struc- 
ture, with two clearly separated peaks for finite B, is not 



yet fully developed. The two peaks largely overlap, and 
the distance of the peaks is below the expected Ae = 2B, 
since tunneling considerably broadens the dot levels. We 
note however, that for larger B, the correct peak spacing 
of Ae = 2B is reproduced (data not shown). 

In order to illustrate the role of the interaction U, wc 
show the dependence of SG on U at eo = in the inset 
of Fig. [TO] For U — 0, both levels contribute e 2 /h to 
the conductance at low temperatures, while for finite U 
these contributions are reduced. The reduction increases 
with growing U, qualitatively consistent with previous re- 
sults. The interaction corrections are smaller away from 
resonance, see inset of Fig. [TOJ but they also grow with 
increasing U. 

Next, we address the temperature dependence of the 
linear conductance (numerically evaluated for eV = 
0.05r). In Fig. EEU we show G(T) for different values 
of U (up to U = AY) at e = B = 0. For U = 1.2Y, the 
deviation from the U = 0-result is small in the consid- 
ered temperature range. For larger U, deviations become 
more pronounced at low temperatures where interaction 
becomes increasingly relevant. Up to present, we have 
obtained converged results in the regime of small bias 
voltages for interaction strengths U < AY for tempera- 
tures above or close to the Kondo temperature, T > Tk- 
The corresponding Kondo temperatures are (see above) 
T K = 0.38r for U = ZY and T K = 0.293r for U = AY. In 
the regime Tk ^ T < 107V, we can compare our results 
to the result of Hamann [60, roll , 



G(T) 



1 



\n{T/T KH ) 



[ln 2 (T/T^) + 37r 2 /4] 1 /2 



(27) 



for the linear conductance, where Tkh = 2V/1.2, see 
Fig. HU (solid lines). In Ref. [60], it has been shown that 
the results of the numerical RG coincide with those of Eq. 
(|27[l in this regime. Fig. [TT] illustrates that the agree- 
ment between the two approaches is satisfactory and 
shows that the ISPI provides reliable results in the linear 
regime above or close to the Kondo temperature. As al- 
ready mentioned, convergence is problematic in the linear 
regime for temperatures lower than the Kondo tempera- 
ture for larger values of U. This implies that the equlib- 
rium Kondo regime is difficult to explore using the ISPI 
approach. However, the situation is more favorable for 
large bias voltages, where short to intermediate memory 
times are sufficient, and we have achieved convergence 
up to U = AY, see Fig. [5] and the next subsection. 



D. Large bias: eV > Y 

Let us then turn to noncquilibrium transport at volt- 
ages eV > r. Here, the transport window is ~ eV, and a 
double-peak structure for dl/dV emerges even for B = 0, 
with distance eV between the peaks. As interaction cor- 
rections to the nonlinear conductance are largest when 
a dot level is in resonance with one of the chemical po- 
tentials of the leads, the double-peak structure is again 
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transferred to S(dI/dV). This can be observed in Fig. 
IT2| where we show results for eV = 3T but otherwise the 
same parameters as in Fig. [TU] For a finite magnetic field, 
each peak of the double-peak structure itself experiences 
an additional Zeeman splitting, resulting in an overall 
four-peak structure. For B = T and the shown values 
of U, the two innermost peaks (closest to eo = 0) over- 
lap so strongly that they effectively form a single peak 
at eo = again. The two outermost peaks are due to 
the combination of the finite magnetic field and the bias 
voltage. 

Increasing the on-dot interaction U leads to a reduc- 
tion of the differential conductance peaks compared to 
the noninteracting case, i.e., the interaction corrections 
are again largest when the level energy matches the chem- 
ical potential in the leads. This is observed in the inset 
of Fig. [12l for B = and eo = ~T, which is close to the 
peak maximum. Away from resonance, the interaction 
corrections are reduced. Note that the four-peak struc- 
ture is already present in the noninteracting case (with 
B =/= 0) and, hence, is not modified qualitatively by U. 

In the regime eV 3> Tk, we can compare our results 
to the perturbative RG result [l2|, HH 

I(V)= 3 -^ln-*(eV/f K ), 

with Tk = 2Tk / s/n. Fig. [5] shows the result for the 
stationary current for U = 4T (where Tk = 0.29r) for 
eV = 2T. The perturbative current amounts to = 
2.28eT/h, while the ISPI value is I IS pi = 2.2heT/h. The 
quite satisfactory agreement suggests that ISPI is indeed 
a reliable new method that holds promise for reaching 
the noncquilibrium Kondo regime. A detailed study of 
the noncquilibrium Kondo effect using ISPI will be given 
elsewhere. 

Finally, we address the temperature dependence of the 
nonlinear conductance dl/dV. ISPI results for eV = 
2T, eo = B = are shown in Fig. [13] Again, as in the 
linear regime, the conductance increases with lower tem- 
peratures, and finally saturates, e.g., at dl/dV = e 2 /h for 
U = and eV = 2T. Clearly the conductance decreases 
when the bias voltage is raised. Increasing U renders this 
suppression yet more pronounced, see also inset of Fig. [13] 
for the corresponding corrections. At high temperatures, 
thermal fluctuations wash out the interaction effects, and 
the interaction corrections die out. 



VII. DISCUSSION AND CONCLUSIONS 

In summary, we have introduced a scheme for the 
iterative summation of real-time path integrals (ISPI), 
and applied it to a prototypical problem of quantum 
transport through an interacting quantum dot coupled 
to metallic leads held at different chemical potentials. 
After integrating out the leads, a time-nonlocal Keldysh 
self-energy arises. Exploiting the exponential decay of 



the time correlations at finite temperature allows to in- 
troduce a memory time t c beyond which the correlations 
can be truncated. Within r c , correlations arc fully taken 
into account in the corresponding path integral for the 
current generating function. Then, through a discrete 
Hubbard-Stratonovich transformation, interactions can 
be transferred to an auxiliary spin field, and an itera- 
tive summation scheme has been constructed to calcu- 
late the transport current. The remaining systematic 
errors due to the finite time discretization and the finite 
memory time t c are then eliminated by a refined Hirsch- 
Fye-type extrapolation scheme, rendering the ISPI nu- 
merically exact. From this construction, it is clear that 
the calculation is reliable for temperatures above a cer- 
tain interaction-dependent temperature scale. The latter 
depends also on the computational power available, and 
vanishes in the absence of interactions. 

The general scheme has been applied to the canoni- 
cal example of an Anderson dot with interaction U, for 
which many features of the transport characteristics are 
well understood. This allows to carefully and system- 
atically check the algorithm. In the regime of linear 
transport, we have recovered results from second-order 
perturbation theory in U in the limit of very small inter- 
action strength, but found significant deviations already 
for small-to-intermediate values of U. In the incoher- 
ent sequential regime, we recover results from a mas- 
ter equation approach. Taking U/T = 4, we have fur- 
thermore reproduced the behavior of the linear conduc- 
tance above the Kondo temperature and found satisfac- 
tory agreement with numerical RG results. In addition, 
we have investigated the regime of correlated nonlinear 
transport, where, in our opinion, the presented method 
is most valuable. We have checked that for large volt- 
age, known results on the noncquilibrium Kondo effect 
are reproduced as well. 

Up to now, we have achieved converged results for 
small-to-intermediate U even at rather low temperatures 
and small bias voltages. However, the strong-coupling 
regime of the Kondo effect has not yet been fully cap- 
tured, as the required memory times become quite large, 
and, at the same time, small 5t are necessary because of 
the large U. This combined requirement makes it difficult 
to reach convergence. Nevertheless, the nonequilibrium 
Kondo regime, representing an intermediate-to- weak cou- 
pling situation, seems tractable by the ISPI scheme. The 
applicability and accuracy of ISPI has been demonstrated 
for bias voltages larger than the Kondo temperature. We 
will present a detailed study of this interesting regime 
elsewhere. 

Our approach is, in fact, similar in spirit to the well- 
established concept of the quasi-adiabatic path inte gral 
(QUAPI) scheme, introduced by Makri and Makarov [57| 
in its iterative version. This method has been developed 
to describe the dynamics of a quantum system coupled 
to a bath of harmonic oscillators held at equilibrium, in 
order to obtain exact results for quantum decoherence 
and dissipation, see also Refs. [H, [62||. For a dipole-type 
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system- bath coupling, as it occurs, e.g., in the spin-boson 
model [HfJ], the bath- induced correlations are encoded in 
the Feynman- Vernon influence functional [uH ]. which is 
solely a functional of a single system operator, say x, 
which couples the system to the bath. The very exis- 
tence of this functional allows to perform a basis rota- 
tion to the eigenbasis of x (called the discrete variable 
representation, DVR). Then, the influence functional can 
be evaluated at the eigenvalues of x during the numer- 
ical calculation of the real-time path integral over x(t). 
Moreover, within the QUAPI scheme, the influence func- 
tional also generates correlations which are non-local in 
time, but which also decay exponentially. This allows to 
truncate them beyond a memory-time r mem , which coin- 
cides conceptually with our correlation time r c . Then, in 
practice, numerical convergence has to be achieved with 
respect to the memory time [53[. Afterwards, the only 
remaining Trotter error can be completely eliminated by 
extrapolation [531 ] . In the present case of noncquilibrium 
transport, several fundamental differences occur, and, in 
fact, only the strategy of memory truncation can be taken 
over. The major difference is that there exists no simple 
Feynman- Vernon influence functional. While the lead 
fermions can be integrated out, see Eq. @, the cor- 
responding Grassmann variables for the dot cannot be 
transferred to a practical computational scheme along 
the lines of the QUAPI method. Instead, here we in- 
tegrate out all fermions, at the expense of introducing 
auxiliary Ising spins via the HS transformation. Then, 
ISPI performs a summation over these Ising spins. An- 
other difference is the form of the tunnel coupling in the 
Hamiltonian. As two non-commuting system operators 



d a and S a occur, no analogue of DVR can be established 
here. 

Finally, we note that we have chosen the Anderson 
model as a simple but non-trivial toy model for quantum 
transport in order to establish and test the numerical 
algorithm. We believe that this approach, or modifica- 
tions thereof, will be of importance in numerical calcu- 
lations of quantum transport properties. Compared to 
other approaches, it has several advantages (e.g., numer- 
ical exactness, direct nonequilibrium formulation, no sign 
problem), but is, on the other hand, also computation- 
ally more costly than most other techniques, especially 
for strong correlations and/or low energy scales (temper- 
ature, voltage). In any case, it goes without saying that 
other interesting models for quantum transport exist, and 
future work will be devoted to apply ISPI to those. 
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FIG. 1: Keldysh contour: every physical time has two con- 
tour representatives on the branches ±. The measurement 
time t m only has a single representative on the upper branch. 





FIG. 2: (Color online) Real part of the Keldysh GF com- 
ponent G~~(t-t') for (a) e = 0,B = 0, T = T, (b) 
e = 0, B = 0, T = O.ir, (c) e = 0, T = O.ir, eV = V, B = 
and B = ST, (d) B = 0,T = O.ir, eV = T,e = T and 
eo = 5r. The inset in (a) shows the absolute values of the 
same data in log-linear representation, highlighting the expo- 
nential decrease. 
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FIG. 3: (Color online) Quadratic dependence of the Trotter 
error as obtained in the convergence procedure for St — > for 
a fixed memory time r c = 0.5/r. Parameters are eo = 0, B — 
0, T — 0.1F. The extrapolated values are given in the figure. 




FIG. 4: (Color online) Dependence of the differential con- 
ductance dI(r c )/dV on the inverse memory time l/r c for 
different values of U after the Trotter error has been elim- 
inated. Parameters are eo = B — 0,T = O.ir. Solid 
fines correspond to a linear fit. The extrapolated values for 
dl/dV — dl(l/r c — > 0)/dV axe given in the figure. 
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FIG. 5: (Color online) Time-dependent current I(t) for two 
representative parameter sets: (1) U — 4F, eV = 2T,T — 
O.ir (black circles) and (2) U = 0.5F, eV = 0.6I\ T = T (red 
squares). Solid line: analytical result for stationary 7 from 
nonequilibrium Kondo theory applied to parameter set (1) 
[13. l5a|. Symbols are ISPI results, dashed lines are guides to 
the eye only, and eo = B = 0. 
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FIG. 6: (Color online) Stationary current as a function of 
eo (a) and of the bias voltage eV (b) for the noninteracting 
case U = 0. Shown are the numerical (circles) and the exact 
analytical (dashed curve) results for (a) eV = 0.2r, and (b) 
e = 0. In both cases, T = 0.1F, B = 0. 
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FIG. 7: (Color online) Interaction corrections SI(V) — 
I(U) — I(U — 0) from the numerical ISPI approach (red sym- 
bols) and from a second-order perturbative calculation (black 
solid lines), for U = O.ir (squares) and U = 0.3F (circles). 
Parameters are eo = B = 0,T = O.ir, and dotted lines are 
guides to the eye only. 



-0.2- 
3, -OA 



-0.6- 



-0.8 





i 1 i 1 


1 1 1 


_ 




""~©— . 
. ~~" G 


1 i 1 i 1 i 1 i 




_ C\J 

«t -0.1 


'Q 




- i°- 2 






to 


1 




0.5 1 1.5 2 

ev/r 

i , i 


i , i 



0.5 



1 

ev/r 



1.5 



FIG. 8: (Color online) Interaction correction SI for the non- 
linear current for U — T, comparing ISPI (red symbols) and 
U 2 perturbation theory (black solid curve). The inset shows 
the corresponding interaction corrections 8{dl /dV) to the dif- 
ferential conductance. Other parameters are as in Fig. [7] 
Dashed lines are guides to the eye only. 



10 



0)_ 

> 

to 



-0.1 



-0.2 



.•0-- 



-0.4 



— master eq. 
O OISPI 



1 I 1 I 


1 1 1 1 1 1 
0- e— e — e — o- 






/ 1 . 1 


. i.i.i 



4 6 

T/r 



8 



4 6 8 10 

T/r 

J i I i L 



10 



FIG. 9: (Color online) Comparison of the ISPI method (red 
symbols) and the master equation approach (solid lines) for 
U = r. Main: Corrections S(dI/dV) of the nonlinear conduc- 
tance as a function of temperature for eV = ST. Inset: Same 
for the linear conductance SG. For the remaining parameters, 
see Fig. [7] 
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FIG. 10: (Color online) Interaction correction SG as a func- 
tion of eo, for U = r, different B, and T = O.ir. Inset: 
Dependence of 8G on U for B — 0, on resonance (eo = 0, red 
circles) and away from resonance (eo = — 1\ black diamonds). 
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FIG. 11: (Color online) Linear conductance G vs temperature 
T, for U = (dashed line) and U = 1.2, 3, 4r (e = B = 0). 
Symbols denote the ISPI results while the solid lines are those 
of Eq. (HZ)). 




FIG. 12: (Color online) Same as Fig. llOl but for the nonlinear 
conductance, taken at eV = 3I\ 
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FIG. 13: (Color online) Nonlinear differential conductance 
dl/dV vs temperature T, for eV = 2F, eo = B = 0, and 
{7 = 0, 1.2, 3r. Inset: corresponding interaction corrections 
S(dI/dV). (The Kondo temperature for U = 3F is T K = 
0.38r.) 



